import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j0, j1
from numpy.fft import fft, ifft, fftfreq
from scipy.linalg import solve

# ====================== 全局配置参数 ======================
config = {
    "a": 0.15,          # 荷载半径 (m)
    "peak": -0.7e6,     # 荷载峰值 (Pa)
    "T_load": 0.6,     # 荷载周期 (s)
    "Nt": 2048,         # 频率采样点数
    "block_size": 1000,   # 分块处理大小
    "dtype": np.complex128,  # 计算精度
    "k_min": 1e-18,     # 波数截断值
    "k_max": 1000.0,    # 波数截断值
    "Nk": 6000,         # 波数采样点数
    "s": 1e6,           # 缩放系数（用于准静态过程，同时缩小弹性模量和密度，提高计算精度，保持量纲稳定）
    "layers": [
        #{'h': 0.18, 'E': 0.4e10, 'nu': 0.25, 'rho': 2300, 'friction': 1},
        #{'h': 0.36, 'E': 0.8e10, 'nu': 0.15, 'rho': 2000, 'friction': 1},
        #{'h': 0.20, 'E': 0.3e10, 'nu': 0.15, 'rho': 2000, 'friction': 1},
        #{'h': 100,  'E': 0.8e8,  'nu': 0.35, 'rho': 1500, 'friction': 0}
        {'h': 0.15, 'E': 0.1e10, 'nu': 0.35, 'rho': 2300, 'friction': 1},
        {'h': 0.25, 'E': 0.2e9, 'nu': 0.35, 'rho': 2000, 'friction': 1},
        {'h': 100,  'E': 0.1e9,  'nu': 0.35, 'rho': 1500, 'friction': 0}
    ]   
}

complex_dtype = np.dtype(config["dtype"])
real_dtype = np.array([], dtype=complex_dtype).real.dtype

# ====================== 材料参数预处理 ======================
def precompute_material_properties(layers):
    """预处理各层材料参数"""
    mask = config["T_load"] >= 1.0
    for layer in layers:
        E = np.where(mask, layer['E']/config["s"], layer['E'])       
        nu = layer['nu']
        rho = np.where(mask, layer['rho']/config["s"], layer['rho'])
        layer['la'] = (nu * E) / ((1 + nu) * (1 - 2 * nu))
        layer['mu'] = E / (2 * (1 + nu))
        layer['cp'] = np.sqrt((layer['la'] + 2 * layer['mu']) / rho)
        layer['cs'] = np.sqrt(layer['mu'] / rho)
    return layers

# ====================== 时间序列生成 ======================
def generate_time_series(config):
    """生成时间序列和荷载频谱"""
    T_total = 61 * config["T_load"]
    Nt = config["Nt"]
    dt = T_total / Nt
    t = np.linspace(0, T_total, Nt)
    
    # 生成荷载信号
    St = 10 * config["T_load"]
    mask = (t >= St) & (t <= St + config["T_load"])
    Ft = np.where(mask, config["peak"] * np.sin(np.pi*(t - St)/config["T_load"]), 0.0)
    
    # 傅里叶变换
    Ft_omega = fft(Ft) * dt
    omega = 2 * np.pi * fftfreq(Nt, dt)
    
    return t, Ft_omega, omega

# ====================== 波数参数计算 (汉克尔变换) ======================
def compute_Q(k, q, a):
    """计算荷载的汉克尔变换"""
    Q = np.zeros_like(k, dtype=real_dtype)
    mask = k != 0
    Q[mask] = (a * q / k[mask]) * j1(a * k[mask])
    Q[~mask] = q * a**2 / 2  # 处理k=0的情况
    return Q

def compute_wavenumber_parameters(config):
    """生成波数并计算荷载的波数分布"""
    k = np.linspace(config["k_min"], config["k_max"], config["Nk"])
    Fr_k = compute_Q(k, 1, config["a"])
    return k, Fr_k

# ====================== 层刚度矩阵计算 ======================
def compute_layer_stiffness(k_block, omega, layer_params):
    """计算单层刚度矩阵"""
    h = layer_params['h']
    mu_val = layer_params['mu']
    la = layer_params['la']
    cp = layer_params['cp']
    cs = layer_params['cs']
    dtype = config["dtype"]
    
    # 维度处理
    k = k_block[:, np.newaxis].astype(dtype)  # 形状 (Nk, 1)
    omg = omega[np.newaxis, :].astype(dtype)  # 形状 (1, Nt)
    Nk, Nt = k_block.size, omega.size

    # 初始化4x4刚度矩阵
    K = np.zeros((Nk, Nt, 4, 4), dtype=dtype)

    # 创建一维布尔掩码 (Nt,)
    omega_zero_mask = np.squeeze(np.abs(omg) < 1e-30)  # 形状 (Nt,)
    non_zero_mask = ~omega_zero_mask
    
    # 初始化分母零掩码 (Nk, Nt)
    denom_zero_mask = np.zeros((len(k_block), len(omega)), dtype=bool)
    
    # ====================== ω=0的特殊处理 ======================
    if np.any(omega_zero_mask):
        λ = la
        μ = mu_val
        kh = k * h  # 广播为 (Nk, 1)
        exp_hk = np.exp(-kh)
        exp_2hk = np.exp(-2*kh)
        exp_3hk = np.exp(-3*kh)
        exp_4hk = np.exp(-4*kh)
        
        # 公共分母
        denom = (
            - (λ+3*μ)**2*(1-exp_2hk)**2
            + 4*(λ+μ)**2*kh**2*exp_2hk
        )

        # K11
        num_K11 = 2*k*μ*(
            - (λ**2 + 5*λ*μ + 6*μ**2)
            + 4*kh*(λ**2 + 3*λ*μ + 2*μ**2)*exp_2hk
            + (λ**2 + 5*λ*μ + 6*μ**2)*exp_4hk
        )        
        K[:, omega_zero_mask, 0, 0] = (num_K11 / denom)

        # K12
        num_K12 = 2*k*μ*(
            - (λ*μ+ 3*μ**2)
            + 4*kh**2*(λ**2+2*λ*μ+μ**2)*exp_2hk
            + 2*(λ*μ+ 3*μ**2)*exp_2hk
            - (λ*μ+3*μ**2)*exp_4hk
        )
        K[:, omega_zero_mask, 0, 1] = (num_K12 / denom)

        # K13
        num_K13 = -4*k*μ*(
            kh*(λ**2 + 3*λ*μ + 2*μ**2)*exp_hk
            - (λ**2 + 5*λ*μ+ 6*μ**2)*exp_hk
            + (kh*(λ**2 + 3*λ*μ + 2*μ**2) + λ**2 + 5*λ*μ + 6*μ**2)*exp_3hk
        ) 
        K[:, omega_zero_mask, 0, 2] = (num_K13 / denom)

        # K14
        num_K14 = -4*h*k**2*μ*(
            (λ**2 + 3*λ*μ + 2*μ**2)*exp_hk
            - (λ**2 + 3*λ*μ + 2*μ**2)*exp_3hk
        )
        K[:, omega_zero_mask, 0, 3] = (num_K14 / denom)

        # K22
        num_K22 = -2*k*μ*(
            + (λ**2 + 5*λ*μ + 6*μ**2)
            + 4*kh*(λ**2 + 3*λ*μ + 2*μ**2)*exp_2hk
            - (λ**2 + 5*λ*μ + 6*μ**2)*exp_4hk
        )
        K[:, omega_zero_mask, 1, 1] = (num_K22 / denom)

        # K24
        num_K24 = 4*k*μ*(
            kh*(λ**2 + 3*λ*μ + 2*μ**2)*exp_hk
            + (λ**2 + 5*λ*μ + 6*μ**2)*exp_hk
            + kh*(λ**2 + 3*λ*μ + 2*μ**2)*exp_3hk
            - (λ**2 + 5*λ*μ + 6*μ**2)*exp_3hk
        )
        K[:, omega_zero_mask, 1, 3] = (num_K24 / denom)
       
        # 对称元素填充
        K[:, omega_zero_mask, 1, 0] = K[:, omega_zero_mask, 0, 1]
        K[:, omega_zero_mask, 1, 2] = -K[:, omega_zero_mask, 0, 3]
        K[:, omega_zero_mask, 2, 0] = K[:, omega_zero_mask, 0, 2]
        K[:, omega_zero_mask, 2, 1] = -K[:, omega_zero_mask, 0, 3]
        K[:, omega_zero_mask, 2, 2] = K[:, omega_zero_mask, 0, 0]
        K[:, omega_zero_mask, 2, 3] = -K[:, omega_zero_mask, 0, 1]
        K[:, omega_zero_mask, 3, 0] = K[:, omega_zero_mask, 0, 3]
        K[:, omega_zero_mask, 3, 1] = K[:, omega_zero_mask, 1, 3]
        K[:, omega_zero_mask, 3, 2] = -K[:, omega_zero_mask, 0, 1]
        K[:, omega_zero_mask, 3, 3] = K[:, omega_zero_mask, 1, 1]               
    
    # ====================== ω≠0的常规计算 ======================
    if np.any(non_zero_mask):
        # 提取非零频率对应的波数和频率
        k_non_zero = k
        omg_non_zero = omg[:, non_zero_mask]
        
        # 计算复数波数（自动广播）
        kp = np.sqrt(k_non_zero**2 - (omg_non_zero**2)/(cp**2) + 0j)
        ks = np.sqrt(k_non_zero**2 - (omg_non_zero**2)/(cs**2) + 0j)
        
        # 预计算指数项（Nk x Nt)
        exp_hkp = np.exp(-h*kp)
        exp_hks = np.exp(-h*ks)
        exp_hkps = np.exp(-h*(kp + ks))
        exp_2hkp = exp_hkp**2
        exp_2hks = exp_hks**2
        exp_2hkps = exp_hkps**2

        # 公共分母
        denominator = (
            (k_non_zero**2 - kp*ks)**2 * (exp_2hkps + 1)
            - (k_non_zero**2 + kp*ks)**2 * (exp_2hkp + exp_2hks)
            + 8 * k_non_zero**2 * kp * ks * exp_hkps
        ) 
        layer_zero_mask = np.abs(denominator) < 1e-300
        denom_zero_mask[:, non_zero_mask] |= layer_zero_mask
        denominator = np.where(layer_zero_mask, 1.0, denominator)  # 避免计算错误（实际会被掩码跳过）
        
        # K11（常规计算）
        numerator = (
            ((-k_non_zero**2 + kp*ks)*exp_2hkps
             + (k_non_zero**2 + kp*ks)*exp_2hkp
             + (-k_non_zero**2 - kp*ks)*exp_2hks
             + k_non_zero**2 - kp*ks)
            * (k_non_zero + ks) * mu_val * kp * (k_non_zero - ks))
        K[:, non_zero_mask, 0, 0] = (numerator / denominator)
    
        # K12
        numerator = (
            -k_non_zero * (-(k_non_zero**2 - kp*ks)*(k_non_zero**2 - 2*kp*ks + ks**2)*exp_2hkps
            + (-12*k_non_zero**2*kp*ks - 4*kp*ks**3)*exp_hkps
            + (k_non_zero**2 + 2*kp*ks + ks**2)*(k_non_zero**2 + kp*ks)*(exp_2hkp + exp_2hks)
            - (k_non_zero**2 - kp*ks)*(k_non_zero**2 - 2*kp*ks + ks**2))
            * mu_val )
        K[:, non_zero_mask, 0, 1] = (numerator / denominator)
    
        # K13
        numerator = (
            2 * mu_val  * kp * (k_non_zero - ks) * (k_non_zero + ks)
            * (exp_hkp*exp_2hks*k_non_zero**2 - exp_2hkp*exp_hks*kp*ks - exp_hkp*k_non_zero**2 + exp_hks*kp*ks))
        K[:, non_zero_mask, 0, 2] = (numerator / denominator)
    
        # K14
        numerator = (
            -2 * mu_val  * k_non_zero * kp * ks * (k_non_zero - ks) * (k_non_zero + ks)
            * (exp_hkp*exp_2hks - exp_2hkp*exp_hks + exp_hkp - exp_hks)) 
        K[:, non_zero_mask, 0, 3] = (numerator / denominator)
    
        # K22
        numerator = (
            -mu_val  * (k_non_zero + ks)
            * ((k_non_zero**2 - kp*ks)*exp_2hkps
                + (k_non_zero**2 + kp*ks)*exp_2hkp
                + (-k_non_zero**2 - kp*ks)*exp_2hks
                - k_non_zero**2 + kp*ks)
            * ks * (k_non_zero - ks))
        K[:, non_zero_mask, 1, 1] = (numerator / denominator)
    
        # K24
        numerator = (
            -2 * mu_val  * ks * (k_non_zero - ks) * (k_non_zero + ks)
            * (exp_hkp*exp_2hks*kp*ks - exp_2hkp*exp_hks*k_non_zero**2 + exp_hks*k_non_zero**2 - exp_hkp*kp*ks))
        K[:, non_zero_mask, 1, 3] = (numerator / denominator)
    
        # 对称元素填充
        K[:, non_zero_mask, 1, 0] = K[:, non_zero_mask, 0, 1]
        K[:, non_zero_mask, 1, 2] = -K[:, non_zero_mask, 0, 3]
        K[:, non_zero_mask, 2, 0] = K[:, non_zero_mask, 0, 2]
        K[:, non_zero_mask, 2, 1] = -K[:, non_zero_mask, 0, 3]
        K[:, non_zero_mask, 2, 2] = K[:, non_zero_mask, 0, 0]
        K[:, non_zero_mask, 2, 3] = -K[:, non_zero_mask, 0, 1]
        K[:, non_zero_mask, 3, 0] = K[:, non_zero_mask, 0, 3]
        K[:, non_zero_mask, 3, 1] = K[:, non_zero_mask, 1, 3]
        K[:, non_zero_mask, 3, 2] = -K[:, non_zero_mask, 0, 1]
        K[:, non_zero_mask, 3, 3] = K[:, non_zero_mask, 1, 1]    
    
    return K, denom_zero_mask

# ====================== 半无限体刚度矩阵 ======================
def compute_halfspace_stiffness(k_block, omega, layer_params):
    """计算半无限体刚度矩阵"""
    mu_val = layer_params['mu']
    la = layer_params['la']
    dtype = config["dtype"]
    
    k = k_block[:, np.newaxis].astype(dtype)
    omg = omega[np.newaxis, :].astype(dtype)
    Nk, Nt = len(k_block), len(omega)

    K0 = np.zeros((Nk, Nt, 2, 2), dtype=dtype)
    
    # 创建一维布尔掩码 (Nt,)
    omega_zero_mask = np.squeeze(np.abs(omg) < 1e-30)  # 形状 (Nt,)
    non_zero_mask = ~omega_zero_mask
    
    # 初始化分母零掩码 (Nk, Nt)
    denom_zero_mask = np.zeros((len(k_block), len(omega)), dtype=bool)

    # ====================== ω=0的特殊处理 ======================
    if np.any(omega_zero_mask):
        λ = la
        μ = mu_val
        denominator = np.array(λ + 3*μ)  # 显式转换为 NumPy 数组后再调整形状
        
        K0_11 = (2*k*λ*μ + 4*k*μ**2) / denominator    # 标量
        K0_12 = 2*μ**2*k / denominator
        
        # 填充到对应位置
        K0[:, omega_zero_mask, 0, 0] = K0_11
        K0[:, omega_zero_mask, 0, 1] = K0_12
        K0[:, omega_zero_mask, 1, 0] = K0_12
        K0[:, omega_zero_mask, 1, 1] = K0_11
        

    # ====================== ω≠0的常规计算 ======================
    if np.any(non_zero_mask):
        cp = layer_params['cp']
        cs = layer_params['cs']
        
        k_non_zero = k
        omg_non_zero = omg[:, non_zero_mask]
    
        kp = np.sqrt(k_non_zero**2 - (omg_non_zero**2)/(cp**2) + 0j)
        ks = np.sqrt(k_non_zero**2 - (omg_non_zero**2)/(cs**2) + 0j)
        denominator = k_non_zero**2 - kp*ks
        
        layer_zero_mask = np.abs(denominator) < 1e-300
        denom_zero_mask[:, non_zero_mask] |= layer_zero_mask
        denominator = np.where(layer_zero_mask, 1.0, denominator)  # 避免计算错误（实际会被掩码跳过）

        # K0[0,0]
        K0[:, non_zero_mask, 0, 0] = ((kp * mu_val * (k_non_zero**2 - ks**2)) / denominator)
        
        # K0[0,1]
        K0[:, non_zero_mask, 0, 1] = ((k_non_zero * mu_val * (k_non_zero**2 - 2*kp*ks + ks**2)) / denominator)
        
        # K0[1,0]
        K0[:, non_zero_mask, 1, 0]  = K0[:, non_zero_mask, 0, 1]  #对称
        
        # K0[1,1]
        K0[:, non_zero_mask, 1, 1] = ((ks * mu_val * (k_non_zero**2 - ks**2)) / denominator)  

           
    return K0, denom_zero_mask

# ====================== 全局刚度矩阵组装 ======================
def assemble_global_stiffness_block(layers, k_block, omega):
    """分块组装全局刚度矩阵"""
    # 确定自由度索引
    layer_indices = [0]
    for i in range(1, len(layers)):
        if layers[i-1]['friction'] == 1:
            layer_indices.append(layer_indices[-1] + 2)
        else:
            layer_indices.append(layer_indices[-1] + 3)
    global_size = layer_indices[-1] + 2
    
    # 初始化分块矩阵
    K_block = np.zeros((len(k_block), len(omega), global_size, global_size), 
                dtype=config["dtype"])
    
    # 初始化全局分母零掩码
    global_denom_zero_mask = np.zeros((len(k_block), len(omega)), dtype=bool)
    
    # 组装各弹性层
    for layer_idx in range(len(layers) - 1):
        layer = layers[layer_idx]
        friction = layer['friction']
        start = layer_indices[layer_idx]
        
        K_layer, layer_denom_zero = compute_layer_stiffness(k_block, omega, layer)
        global_denom_zero_mask |= layer_denom_zero
        
        if friction == 1:
            # 完全粘结条件：填充4x4上三角
            for i in range(4):
                for j in range(i, 4):
                    K_block[:, :, start+i, start+j] += K_layer[:, :, i, j]
        else:
            # 滑动条件：特殊填充模式
            # 填充3x3上三角
            for i in range(3):
                for j in range(i, 3):
                    K_block[:, :, start+i, start+j] += K_layer[:, :, i, j]
                # 填充前三行的第四列
                K_block[:, :, start+i, start+3] += K_layer[:, :, i, 3]
            # 填充第四行第四列
            K_block[:, :, start+3, start+3] += K_layer[:, :, 3, 3]
    
    # 组装半无限体
    last_layer = layers[-1]
    K_halfspace, half_denom_zero = compute_halfspace_stiffness(k_block, omega, last_layer)
    global_denom_zero_mask |= half_denom_zero
    start = layer_indices[-1]
    K_block[:, :, start:, start:] += K_halfspace
        
    # 填充下三角
    for i in range(global_size):
        for j in range(i+1, global_size):
            K_block[:, :, j, i] = K_block[:, :, i, j]
            
    return K_block, global_denom_zero_mask

# ====================== 分块求解器 ======================
def solve_block(K_block, Fr_block, Ft_omega, omega, denom_zero_mask):
    """分块求解位移响应"""
    block_size, Nt = K_block.shape[0], K_block.shape[1]
    u_block = np.zeros((block_size, Nt), dtype=config["dtype"])
    omega = omega.astype(real_dtype)
    
    for m in range(block_size):
        for n in range(Nt):
            # 跳过分母为零的情况
            if denom_zero_mask[m, n]:
                u_block[m, n] = 0.0
                continue
                
            try:
                # 构建荷载向量（作用在第二个自由度）
                rhs = np.zeros(K_block.shape[-1], dtype=config["dtype"])
                rhs[1] = Fr_block[m] * Ft_omega[n]
                
                # 求解线性系统
                #sol = solve(K_block[m, n], rhs)  #精度不高
                sol = np.linalg.lstsq(K_block[m, n], rhs, rcond=1e-300)[0]
                u_block[m, n] = sol[1]
            except np.linalg.LinAlgError:
                u_block[m, n] = 0.0
                
    return u_block

# ====================== 主计算流程 ======================
def main():
    config["layers"] = precompute_material_properties(config["layers"])
    t, Ft_omega, omega = generate_time_series(config)
    Nt = len(t)
    dt = t[1] - t[0]
    
    # 计算波数参数
    k, Fr_k = compute_wavenumber_parameters(config)
    Nk = len(k)  # 确保Nk与实际一致
    dk = k[1] - k[0]  # 波数间隔

    # 初始化完整的频域响应数组 (Nk, Nt)
    w_k_omega = np.zeros((Nk, Nt), dtype=config["dtype"])  #(Nk, Nt)    
        
    # 分块处理主循环
    block_size = config["block_size"]
    total_blocks = (Nk + block_size - 1) // block_size
    
    for block_idx in range(total_blocks):
        start = block_idx * block_size
        end = min(start + block_size, Nk)
        k_block = k[start:end]
        Fr_block = Fr_k[start:end]
        
        print(f"Processing block {block_idx+1}/{total_blocks} (k={start}-{end-1})")
        
        # 分块组装刚度矩阵并求解
        K_block, denom_zero_mask = assemble_global_stiffness_block(config["layers"], k_block, omega)
        u_block = solve_block(K_block, Fr_block, Ft_omega, omega, denom_zero_mask)

        # 将当前块的u_block填充到w_k_omega中
        w_k_omega[start:end, :] = u_block        

        # 释放内存
        del K_block, u_block


    # 预计算贝塞尔函数矩阵 J0(kr) (形状: Nr, Nk)
    r_points = np.array([0, 0.2, 0.3, 0.45, 0.6, 0.9, 1.2, 1.5, 1.8])
    Nr = len(r_points)
    J0_kr = j0(r_points[:, np.newaxis] * k[np.newaxis, :])  # (Nr, Nk) 
    w_r_omega = np.zeros((Nr, Nt), dtype=config["dtype"])  # (Nr, Nk) 

    # 数值积分逆汉克尔变换
    for i in range(Nr):
        # 轴向位移分量 (J0积分)
        integrand_w = J0_kr[i, :, np.newaxis] * w_k_omega * k[:, np.newaxis]
        w_r_omega[i, :] = np.trapz(integrand_w, x=k, axis=0)    # 沿波数维度求和    
    
    # 逆傅里叶变换得到时域结果
    w_r_t = np.real(ifft(w_r_omega / dt, axis=1)) 
    mask = config["T_load"] >= 1.0
    w_r_t = np.where(mask, w_r_t/config["s"], w_r_t)
    
    # ====================== 结果可视化 ======================
    # 创建画布（增加高度到8英寸）
    fig = plt.figure(figsize=(6.5, 5), dpi=300)  # 高度从7增加到8英寸
    colors = plt.cm.jet(np.linspace(0, 1, len(r_points)))
    
    # ========================= 时程图 =========================
    ax1 = plt.subplot2grid((5, 1), (0, 0), rowspan=3)
    
    # 绘图设置
    for i, r_val in enumerate(r_points):
        ax1.plot(t, w_r_t[i], color=colors[i], label=f'r={r_val:.2f}m', lw=1)
    ax1.set_title(f'Hankel Transform($T_d$={config["T_load"]}s)', pad=12)
    ax1.set_xlim(9.5*config["T_load"], 12.5*config["T_load"])
    ax1.grid(True, alpha=0.3)
    ax1.set_ylabel('Displacement (m)', labelpad=-65)
    ax1.set_xlabel('Time (s)', labelpad=-30)
    
    # ========================= 弯沉盆图 =========================
    ax2 = plt.subplot2grid((5, 1), (3, 0), rowspan=2)
    
    # 计算最小弯沉值（最大沉降）
    min_w = np.min(w_r_t, axis=1)
    
    # 绘图设置
    ax2.plot(r_points, min_w, 's-', color='maroon', markersize=4)
    for i, r_val in enumerate(r_points):
        w_val = min_w[i]
        ax2.text(r_val - 0.03, w_val - 0.0001, f'{w_val*1000:.2f}mm', fontsize=10)
    
    ax2.set_ylim(np.min(min_w)*1.3, np.max(min_w)*0.3)  # 自动调整Y轴范围
    ax2.grid(True, alpha=0.3)
    ax2.set_ylabel('Max Deflection (m)', labelpad=-70)
    ax2.set_xlabel('Radial Position (m)', labelpad=-30)
    
    # ======================= 图例设置 ========================
    handles, labels = ax1.get_legend_handles_labels()
    fig.legend(handles, labels,
              loc='lower center',
              ncol=5,
              borderaxespad=0.5,
              columnspacing=1.5,
              bbox_to_anchor=(0.5, 0.03),
              frameon=False)
    
    # ======================= 布局调整 ========================
    plt.subplots_adjust(
        hspace=0.4,
        bottom=0.18,
        top=0.92
    )
    plt.savefig(f'layered_displacement_response_1(T={config["T_load"]}s).png', bbox_inches='tight')
    plt.show()

if __name__ == "__main__":
    main()